Ground state correlations in a trapped quasi one-dimensional Bose Gas 
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We review the basic concepts of a non-equilibrium kinetic theory of a trapped bosonic gas. By extending 
the successful mean-field concept of the Gross-Pitaevskii equation with the effects of non-local, two particle 
quantum correlations, one obtains a renormalized binary interaction and allows for the dynamic establishment 
of non-classical many-particle quantum correlations. These concepts are illustrated by self-consistent numerical 
calculations of the first and second order ground state quantum correlations of a harmonically trapped, quasi 
one-dimensional bosonic gas. We do find a strong suppression of the density fluctuations or, in other words, an 
enhanced number squeezing with decreasing particle density. 
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I. INTRODUCTION 

Already the state of a single quantum mechanical particle 
can exhibit an amazingly complexity. Without the restriction 
to one space dimension or by exploiting additional symme- 
tries, it is virtually impossible even to visualize the informa- 
tion that is encoded in a single particle wave function. Even 
more so, this situation gets quickly out-off hand by adding 
more interacting particles to the system. While the new field 
of quantum information science is exactly trying to accom- 
plish this task of describing and preparing many-particle en- 
tanglement, it remains a formidable challenge. In contrast, 
many-particle physics is blessed with serendipity as the astro- 
nomic growth in the dimensionality of the state space leads 
to a great reduction of complexity. No macroscopic physi- 
cal phenomenon does rely on a miniscule detail of a particular 
quantum state as long as the whole ensemble exhibits a certain 
macroscopic characteristic. 

The past decade of research in degenerate atomic gases has 
produced an amazing wealth of condensed many -particle phe- 
nomena [1^]: Bose-Einstein condensation [5-7], the creation 
of vortices [8, 9], Abrikosov lattices [10], the Mott phase tran- 
sition [11-13], the creation of molecular condensates [14, 15], 
the BCS-BEC crossover [16, 17] and the one-dimensional 
Tonks-Girardeau gas [18, 19], most of which have been seen 
or, at least, were predicted in other, traditional fields of low 
temperature physics. However, observing the original BEC 
phase transition of a dilute atomic gas in situ and real space is 
still a delightful lesson in fundamental physics, as one can be 
an eye witness of the establishment of off-diagonal long range 
order (ODLRO) [20, 21] and get an impression of role that is 
carried by quantum fluctuations . 

The art of performing successful many-particle calcula- 
tions consists of picking the right approximation scheme that 
matches the experimental system on one hand and that is the- 
oretically tractable on the other hand. In the context of weakly 
correlated dilute Bose gases, the Gross-Pitaevskii (GP) mean- 
field picture has been a tremendously rewarding concept and 
the extension to incorporate quantum fluctuations dynami- 
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cally is in principle straight forward, although it involves in 
detail some intricate calculations [22-29]. The description of 
the temporal relaxation [30, 31], the buildup of spatial corre- 
lation functions [32, 33] and the squeezing of atomic number 
density fluctuations [34] follow, consequently. 

However, it has been recognized very early on that the spa- 
tial dimensionality of a system is of utmost importance to its 
physical behavior In particular, it has been proven that the re- 
duction of the available phase space volume leads to enhanced 
fluctuations and the absence of ODLRO in one and two di- 
mensions [35, 36]. During the last years, this fascinating ob- 
servation has received much attention [37^0] as trapped, in- 
homogeneous systems violate the translational invariance as 
required in [35, 36]. This different physical response should 
also be reflected by tuning the trap geometry dynamically 
from 3d^2d^ quasi Id Id [18, 19, 41^3]. The extended 
mean-field theory that will be presented in the following, is 
well suited to describe the cross-over physics with the excep- 
tion of the strongly interacting Tonks-Girardeau regime. 

We have arranged this article according to the following 
outline: In Sec. II, we will give a basic review of the premises 
and concepts of non-equilibrium kinetic theory. In particular, 
we discuss the physical meaning of the relevant master vari- 
ables that are used. We discuss the derivation and approxima- 
tions that lead to the the basic self-consistent time-dependent 
Hartree-Fock-Bogoliubov (HFB) equations of motion in the 
absence of collision. A short review on the mathematical 
properties of "Bogoliubov-like" symplectic self energy oper- 
ators is given, taking special care of the presence or absence 
of a zero energy mode. Sec. Ill is devoted to an application 
of the general formulation to a quasi one-dimensional trapped 
bosonic gas. To calculate specific numbers, we are assum- 
ing in here the typical data of a ^^Rb experiment. In partic- 
ular, we calculate for the zero-temperature ground state of a 
gas: the mean-amplitude, the quantum depletion, the pairing 
field, as well as the first and second order correlation func- 
tions for a full range of particle numbers N ~ (10°, . . . , 10^). 
In Sec. IV, we do draw conclusion and give an outlook to the 
work in progress. Finally, five short appendices compile some 
technical methods or basic statements that were used in the 
article. 
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II. COLLISIONLESS KINETIC EQUATIONS 
A. Quantum dynamics 

The kinetic evolution of a trapped atomic gas is described 
very well by a dilute gas Hamiltonian [1]. In the Umit of 
strongly rarefied atomic gases, it consists primarily of the sin- 
gle particle energy of atoms in a harmonic trap potential Ho 
and the mutual interaction energy amongst all pairs of atoms 
which is mediated through a short-range inter-atomic poten- 
tial Vbin- Strong collisions between atomic triples are very 
unUkely events in this dilute gas limit and can be disregarded, 
consequently. Thus, one finds for the Hamilton operator 

H = J (fx al (^-^A + Ho(x) - 

+^ / '^^^y ^i^y ^'"'"(^ ~ ^y^'" 

where m denotes the atomic mass and is a conveniently cho- 
sen zero-energy reference that will be identified later with the 
chemical potential. 

In the language of second quantization, the action of a field 
operator or on a state in Fock space represents the re- 
moval or creation of another unstructured particle at the spa- 
tial position X. The bosonic nature of the indistinguishable 
particles is reflected by the commutation relation of the fields 
as 



(2) 



Due to the symmetry [a^^jay] = of the bosonic particles 
under coordinate exchange, only even parity contribution of 
the interaction potential Vbin(x) — Vbin(— x) contribute to the 
kinematic evolution. 

In principle, all dynamic and static aspects of the evolu- 
tion of observables O can be obtained from the solution of 
Heisenberg's equation dO/dt = i [H, 0]/h and the knowl- 
edge of the initial state of the system, which is represented 
by the many-body density matrix p. As all observables are 
formed by the elementary quantum fields a , it is only neces- 
sary to consider the Heisenberg equation 

ih^a^{t) = W(x) «x + y d^y Hin(x - y) aJ^Oya^. (3) 

For convenience and later reference, we have introduced in 
here the Hamilton operator of a single trapped atom 



W(x) = -|-A + yho(x)-M. 



(4) 



In principle, any complete state representation can be used 
to perform further calculations. However, it is intuitively clear 
that a representation that matches the geometry better or in- 
corporates conserved symmetries of the system will lead to 
a greatly simplified description, reveal the essential physics 
more clearly and make numerical simulations efficient. Thus, 



we will decompose the quantum field in the position represen- 
tation |x) 

in another complete, yet unspecified basis { 1 1) s \ii,ji,ki)}, 
which is supposedly more suitable. Furthermore, we will also 
employ an implicit summation convention over the quantum 
labels necessary to specify a state completely. This means that 
summation symbols are omitted and a repeated occurrence of 
a dummy summation index on one side of an equation impUes 
a summation. 

In this generic basis, the dilute gas Hamiltonian of Eq. (1) 
reads 



H = Ti}'^ ala^ + (j)^'^^'^ a\ala^a^, (6) 
where we have introduced the matrix elements 



^ {l\n\2) = J d3x(l|x)H(x)(x|2), (7) 

of the single particle Hamiltonian Eq. (3) and the two-particle 
matrix elements 



il234 



-(5)(l|«)(2|Fbm(xi-X2)|3)«)|4) 
^1243 ^ ^2134 ^ ^2143 

2 3 



1 



(8) 



from the binary interaction potential. Due to the bosonic na- 
ture of the particles, only the symmetric part (S) of the in- 
stantaneous coupling vertex cp^'^^'^ is physically relevant and 
the diagrammatic representation also carries this property. 

In the low kinetic energy range that we are interested in, 
repulsive s-wave scattering is the dominant two-particle scat- 
tering event [44, 45]. Provided that a proper T-matrix scat- 
tering calculation has been performed [39, 46] or one has ob- 
tained the experimental scattering data, one can encode this 
information efficiently via a pseudo potential method [47, 48]. 
In the most elementary invocation of the method, one uses 
a fictitious contact potential T4iin(xi,x2) — VoS{xi — X2) 
with a single parameter Vq. This parameter is directly re- 
lated to the scattering length Us of two particles in vacuo by 
Vq = Anfi^as/m, provided we would Umit all physical ap- 
proximations to a first order contribution in Vo- Otherwise 
again an infinite order resummation takes place and will lead 
to a renormalization of the effective coupling constant Vq. We 
will demonstrate this renormalization of the effective coupling 
constant takes place in a self-consistent calculation and also 
leads to a natural momentum cut-off. In the case of such a 
contact potential, one finds for the two-body matrix elements: 



/1234 



^ I d=^ar(l|x)(2|x)(x|3)(x|4), (9) 
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which need not be symmetrized, as it is symmetric already. 

In general, no exact solutions of the field equation Eq. (3) 
or eigen-states of the Hamiltonian Eq. (1,6) are known and the 
few celebrated exceptions [49] such as the Tonks-Girardeau 
gas [18, 50], the Lieb-Lininger solution [37], the Richardson 
pairing model [51, 52] and fermionic Luttinger liquids [53] 
serve as testing grounds to prove the approximation schemes. 
Fortunately however, most physical phenomena are of univer- 
sal character. Thus the system under investigation does not 
have to follow precisely a particular model to show a certain 
response and various approximations are admissible as long 
as the main universal aspects of the problem are incorporated 
in the specific model hypothesis. 



B. Reduced state description with master variables 

The method of quasi-averages and self-consistent field 
equations has a long standing tradition in the description of 
classical gases and fluids [54], in plasma physics, nuclear mat- 
ter physics [55] and condensed matter physics [56-59]. It has 
been applied successfully to classical particles and degenerate 
superfluid bosonic as well as fermionic systems. 

The basic premises for a reduced state description of a 
weakly correlated many-body problem relies on the existence 
of a well separated hierarchy of time, energy and length 
scales. If this is the case, one can assume that the information 
required to describe an ensemble effectively can be parame- 
terized with a set of a few relevant variables [56]. However 
the price that has to be payed for reducing the "astronomical" 
dimension of the linear many-body Schrodinger equation is 
giving up the superposition principle and embarking on non- 
linear mathematics. 

In the mean-field approximation, we want to assume that 
the ensemble of relevant quantum states establishes a well de- 
fined mean value for the field operator and that quantum fluc- 
tuations only cause small deviations around it 

«! = {ai)+Sa^, (10) 

such that {6di) = 0. As mentioned before in Eq. (5), we use 
in here the shorthand notation for any complete set quantum 
of labels e.g., 1 = xi in position space or 1 = fci in momen- 
tum space, respectively. This number symmetry breaking ap- 
proximation is tremendously useful and can be envisaged also 
as the semi-classical limit of coherent many-particle quantum 
physics in analogy to the description of the optical laser [60]. 
However, one must be aware of the implied consequences that 
have been discussed in the Uterature [2, 23, 26, 61] and num- 
ber conserving approximations [62-64] have their merits, but 
shortcomings as well. 

Based on these assumptions, we want to sum the c-number 
field amplitudes = (a^) that were introduced in Eq. (10) 
over this complete set of states and form a basis-independent 
element of a Hilbert space as 

|a) = ai |1) = (a) = ™- . (11) 



In here, we have used the Dirac notation with the implied 
standard scalar product (a|/3) = / d^x q;*(x)/3(x). Due to 
the non-linear nature of the ensuing mean field equations, we 
can no longer rely on the dynamical superposition principle 
available in linear quantum mechanics. However, the advan- 
tage of the notation arises from preserving and emphasizing 
the geometrical transformations properties of aU correlations 
functions under a change of basis or, most generally, frame of 
reference. The wiggly line has been introduced to represents 
the mean-field amplitude graphically. This symbol hterally 
denotes the state at time t and not the time-ordered propaga- 
tor that evolves it. 

In an analogous fashion, we can separate the single particle 
density operator of the atontic gas / = (a^a) = f^''^ + f 
into a mean-field contribution f^'^'^ and a fluctuation / around 
it. This contribution / is also known as the normal quantum 
depletion of the atomic cloud. Both quantities are hemutian 
tensor operators of rank (1,1) and defined as 

/(^) = |a)(a|, / = /i4 |1)(4|= (12) 

For example, in a position representation, the single parti- 
cle density matrix reads as (xi|/|x2) = /^'^^(xi,X2) -|- 
/(xi,X2) = q:(xi)q!*(x2) + {{al^ - a*{x2))(a^^ - a{xi))) 
Similarly, we define anomalous averages or pairing fields 
m = {aa) = m^"^ + rh, as synmietric tensors of rank (2,0), 
i.e. , 

m^^'^ = \a) \a) , m = mi2 |1) |2) = -h- , (13) 

as weU as their symmetric conjugates as to* = > « . 
These anomalous averages are also associated with the quan- 
tum depletion of the degenerate gas. In contrast to the sin- 
gle particle nature of normal density matrix of Eq. (12), it 
can be seen immediately from the structure of Eq. (13) that 
the pairing field is a two-particle state. It will be shown in 
Sec. II F that the pairing field to(xi,x2) evolves basically 
Uke a generalization the bare two-particle Schrodinger state. 
Thus, it carries all the important physics of binary scatter- 
ing. In a position representation, the pairing field is given by 
(xi| (X2I m = m('=)(xi,X2) +to(xi,X2) = a(xi)a(x2) -|- 
((ax2 - a(x2))(axi - 

C. Structure of the generalized density matrix 

From the transformation properties of the matrices under 
coordinate change, one finds that the normal density ma- 
trix f{t) and the pairing fields m{t) are not independent but 
actually the components of a generalized density matrix G 
[55, 58, 65]. If we represent the system in a state space of 
dimension n, then it is convenient to arrange the mean field in 
a 2n-dimensional row vector x and the fluctuations as a posi- 
tive semi-definite matrix dim[G'] — 2n x 2n, 
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The non-negativity of this co-variance matrix implies that the 
magnitude of the anomalous fluctuations is limited by the nor- 
mal depletion through a Cauchy-Schwartz inequality (see Ap- 
pendix A). In the general context of Green function's, this 
single-time density operator G{t) can also be viewed as a par- 
ticular limit of a time-ordered (T), two-time Green function 
G(T,t), i.e., G{t) = G>{t) = limr^t+T G{T,t). Conse- 
quendy, it is also necessary to consider the opposite limit and 
to define a time-reversed, single-time density operator through 
G<(t) = lim^^4_ TG{T,t). Explicitly, this operator is given 
by 



= aiG>Vi = G^ -h(73 = 



1 + / m 
m* f* 



(15) 



where standard Pauli spin matrices have been introduced and 
are defined in Appendix B. 

The specific structure of the generalized density matrix im- 
plies various important physical properties. First of all, we 
have to assume that there is a basis that diagonalizes this fluc- 
tuation matrix. Exactly n of its 2 n eigen- values correspond to 
the positive occupation numbers of finding a particle or, more 
generally, a quasi-particle in a certain mode. For a given, but 
otherwise arbitrary, G matrix, one can construct this basis by 
studying the transformation law of the density matrix under a 
canonical transformation T (see Appendix C), 



G' = TGTK 



(16) 



It is important to note that this is not the transformation law of 
a general matrix under coordinate change. This would require 
that = T~^. However, by only using the properties of 
the symplectic transformations, one can show that a canonical 
eigen- value problem is defined by 



(0-3 G) Tt = (era G') 



(17) 



The solution of this eigen-value problem yields the eigen- 
vector matrix and the corresponding diagonal eigen-value 
matrix as G'. All normalizable states can be rescaled such that 
T(jz =0-3. Now, we are able to reconstruct the positive G 
matrix 



G = VPV\ 



(18) 



from its eigen-vectors V = a^T^ and the diagonal, positive 
occupation number matrix P = as G' as- 
Second, an important feature of an admissible fluctuation 
matrix is its consistency with the conomutation relation, i. e. , 

5i2 and (0.1(12) = (S2^i)- This was al- 



''l"-2/ 



ready expressed in Eq. (15) as 

criG*(j\ — G = (73. 



(19) 



By invoking the properties of a unitary symplectic trans- 
formation, one can show that the elements of the diagonal 
occupation number matrix P are not 2n independent vari- 
ables. Actually half of them are determined by the other half, 

P{n+l,...,2n) = 1 + P(l,...,n)^ Or 

aiPai-P = (73. (20) 



In other words, by separating the occupation numbers P and 
the eigen-vector matrix V into a first and second half, i. e. , 
P = {P+,l+ P+) and V = {V+,V-), one can then decom- 
pose a general fluctuation matrix as 



G =V+P+Vl + V_ (1 + P+) Vl. 



D. Dynamic equations of motion 



(21) 



It is now straight forward to derive equations of motion for 
these averages directly from the Heisenberg equation Eq. (3). 
However, due to the nonUnear character of the operator equa- 
tion, one finds always a coupling of correlation functions in- 
volving n fields to a correlation function of n + 2 fields. 
If we truncate this infinite hierarchy of correlation func- 
tions (BBGKY) at the level of one and two-operator fields 
{a , /, m} and approximate higher correlation functions with 
the help of Wick's theorem [27, 28, 56], one obtains the fol- 
lowing equations of motion, also known as time-dependent 
Hartree-Fock Bogoliubov (THFB) equations (THFB) 



ih^G 

dt 



= nx, 

= SG-GSt. 



(22) 



(23) 



For the evolution of the mean-field oris finds a general- 
ized Gross-Pitaevskii propagator that is defined as 



n 



Ha/- 



-n* 



-n 



(24) 



The two contributions that define this symplectic propagator 
are a normal hermitian Hamiltonian operator and a symmetric 
anomalous couphng potential 



— + I +2 X 



(25) 



o 



(26) 



It is easy to identify with the well known hermitian GP- 
propagator that accounts for the free evolution of the mean- 
field Ti, its self-interaction U^ic), as well as the energy shift 
Up which caused by the presence of the non-condensate 
cloud. However, due to the existence of the anomalous fluctu- 
ations there is also a coupling through 11^ to the time-reversed 
field. For convenience, we have introduced two auxiliary po- 
tentials Uj and V„, . Explicitly, they are defined in terms of 
the two-body matrix elements as 

Zif = 2(/.1234/32, (27) 

V^2=2</,i23^m34. (28) 
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In a position representation, this reduces to the familiar non- 
local Hartree-Fock potentials 

^/(x, y) =^ [Hin(x - y) /(x, y) (29) 

+^(x - y) y dh Fbin(z - y) /(z, z)], 

Vrn (x, y ) =ybm (x - y ) TO (x, y ) . (30) 

Similarly, one finds that the evolution of the density oper- 
ator G is ruled by a HFB self-energy E that can be obtained 
also by variational methods [58]. In detail, this symplectic 
self-energy is given by 

S = ( ^ ) , (31) 

where we have introduced hermitian Hamiltonian operators 
and symmetric anomalous coupling potentials as 

Sa/- = H + 2Uf(c)+2Uf (32) 




It is important to note the different weighing factors of the 
mean-field potential in Eqs. (25) and (32) and the fact that the 
potentials are local in time. 



the indefinite norm = W^a^^W, except for those that 

belong to zero or complex eigen- values. It is important to 
note that this energy basis W is in general distinct from the 
instantaneous basis V that diagonaUzes the fluctuation matrix 
G in Eq. (18). They do coincide only in equilibrium. The 
mathematical properties of the eigen-states W can be derived 
easily from the intrinsic symmetries of the HFB self-energy 
operator: 

S = -cTiSVi, (35) 
Et = asEag. (36) 

Thus, if is a solution of the right eigen-value problem with 
energy E, 

i:W = WE, (37) 

it follows directly from Eq. (35) that W = ai W*, is also 
a right eigen- vector but corresponds to the eigen-value E = 
—E*. Starting from the second symmetry in Eq. (36) and 
the right eigen-value problem of Eq. (37), it is easy to con- 
struct the left eigen-vectors W = that correspond to 
the eigen- values E = E*: 

Wi: = E*W. (38) 

Finally, from a combination of the results for the right and left 
eigen-vectors, it follows that the eigen-vectors are orthogonal 
with respect to the metric : 

Q = {E* -E')WlasWE,, (39) 

if E* ^ E'. On the other hand, this relation implies also that 
eigen-vectors that belong to complex eigen- values must have 
zero norm. 



E. Structure of the Hartree-Fock-Bogoliubov self energy 

1. Normal quasi-particle modes 

Symplectic self-energy operators arises not only naturally 

in kinetic theories [28, 66] or variational calculations, but in 
many other contexts involving stability analysis. In the case of 
bosonic fields, the self-energy operator is of the generic form: 

In here, E ^ stands for a hermitian operator E^^^ = Ej^ and 
E^ denotes an anomalous coupling term that has to be sym- 
metric E^ = Ey^. The relative size of the operators Eyv" 
and E^ determines the character of the energy spectrum. It 
can either be real valued with pairs of positive and negative 
eigen-energies, or one finds a doubly degenerate zero eigen- 
value, if the energy difference between the smallest positive 
and highest negative vanishes (gap-less spectrum). In the gen- 
eral case, there is a mixed spectrum consisting of pairs of real 
sign-reversed as well as pairs of complex conjugated eigen- 
values. The eigen-vectors W are normalizable with respect to 



2. Defective sub-spaces 

The situation of a doubly degenerate zero energy eigen- 
value E = (] needs special attention. One can view this case as 
a limit when two non-degenerate states approach each other. 
However, as the energy gap decreases, the two eigen-states 
become more and more collinear. Thus, in the limit of a van- 
ishing energy separation, the dimension of the spanned vector 
space collapses from 2 to 1 and E becomes defective [58, 67]. 
In the present context, a gap-less linear response matrix occurs 
from a perturbation analysis of the simple Gross-Pitaevskii 
equation and describes the collective excitation of the system. 
The emerging zero energy state is called Goldstone mode and 
can be interpreted physically as an attempt to restore the bro- 
ken number symmetry. On the other hand, we find for the 
THFB Eqs. (22,23) that the self-energy matrix E has a gap 
and its eigen-states form a complete non-defective basis. It is 
important to distinguish these states from the collective exci- 
tations of the total system, whose excitation spectrum is again 
gap-less [29, 68]. 

The general situation can be described by separating the 
two ground state modes {WEQ,W-Ea) from the remaining 



6 



states W such that dim[W^']=2n x 2(n — 1) and the diagonal 
eigenvalue matrix has a dim[i<^']=2(n — 1) x 2 (n — 1). By in- 
troducing two quadrature modes PQ with dmi[PQ]=2n x 2 



via 



SW = W'E', 
= PQih 







(40) 
(41) 



one can span a two dimensional vector space which is orthog- 
onal to the higher modes W and does not collapse [69, 70]. 
This is mathematically achieved by the construction of a 
"best" basis in the context a singular value decomposition [67] 
and introduces two singular values, i. e. , an inertial mass pa- 
rameter M and a gap energy Eq = hu>o of the number and 
phase quadratures. The quadrature states satisfy the following 
orthogonality 



PQV3VK' = 0, PQ^a^PQ ^ ha. 



(1) 



(42) 



and symmetry relations PQ = —aiPQ*. From a dimen- 
sional consideration, it is obvious that one also has to use 
lower dimensional Pauli matrices, which act in the appropriate 

subspaces such that dim[(7[2,i 2 31=2/ x 21. All states together 
form again a complete basis such that 



(1) 



W'ai^'-'^W'^ +PQ\- PQ^ = as. 



h 



(43) 



With these states, we can then obtain the following spectral 
decomposition of the self-energy 



Eas = W'{E'at'^)W'' + PQ " j PQt 

The physical meaning of the quadrature states can be under- 
stood most clearly when mapping them again back onto the 
quantum field 







(45) 



Now, these field quadratures Q, P satisfy a conventional 
Heisenberg-Weyl algebra and define in depend, as well as or- 
thogonal bosonic quasi-particles 6j according to 



[Q,p]=ih, 



(46) 



F. Upgrading off-diagonal potentials to many-body T-matrices 

While one can understand the physical structure of the 
self-energy best in the general form of the kinetic equations 

Eq. (23), one can appreciate other aspects much better by con- 
sidering the equations for the components individually. 



dV 



mil 



(47) 



in—m 
dt 



m + m E> + (1 + /)* + /S^, (48) 



In stationarity and by assuming a real valued self-energy, one 
can solve for m in an HF eigen basis 



and finds 



2/i - (ei -I- £2) 



mi2 



(49) 



(50) 



No rules for circumventing poles in the above energy denom- 
inator need to be specified as all energies satisfy ei> ji. 

In the present article, we have only considered first order 
processes and found the one-loop contribution to the self- 
energy as 



,(<=) 



+ 



A) 



(51) 



Now, if we substitute the stationary solution of Eq. (50) into 
this equation, one obtains an inhomogeneous linear equation 
for the self-energy 



1,1234 



(]l + 2/)33' , 

2/i - (e3 + ti) ' 



^3'4 



(52) 



It is now clear that the value of T,j^ is only determined by 771*^°^ 
in the absence of a homogeneous contribution. Thus, there 
must be a T-matrix that maps mean-field to the off-diagonal 
self-energy 



t12 



t1234/ 



34 



(53) 



1 



This leads to the following general definition for an off-the 
energy shell T-matrix 

T^^\e) = 2(/.1234 + 2<^1^3'4' (l + 2/)3'3" 3"4'34_ ^54) 

e — (es' + £4') 
In operator notation this is equivalent to 

TA{e) = V + VGA{e)[{t + 2/) ® t\Tj^{e). (55) 

where Gj^{e) = (e — E^/- (g) 1 — 1 E^/-)^-' is a Green's 
function, which describes the propagation of two independent 
HF particles according to Eq. (49). The diagrammatic repre- 
sentation of this T4-matrix relation is 




(56) 



7 



G. Invariants and Conservation laws 

1. Number 

The total particle number N = J (Px at.a^ is a conserved 
quantity if the atoms evolve under the generic two-particle 
Hamiltonian operator H given by Eq. (1), i. e. , [H, N] = 0. 
This conservation law implies that the system is invariant un- 
der a global phase change d ^ a exp (— i $). By using this 
continuous symmetry, i. e. , a a exp (— i $), / /, and 
fh ^ fh exp (— 2 i $), it is easy to see that the kinetic equa- 
tions Eqs. (22,23) are also explicitly number conserving at all 
times: 

N{t) = (N) = Tr + m} = const. (57) 

Nevertheless, it is important to note that out-of equilibrium 
there can be a continuous particle exchange between the con- 
densate and the non-condensate clouds. 



2. Energy 

In the absence of any time-dependent external driving 
fields, such as optical lasers or magnetic rf-fields, the overall 
energy H must be conserved as well. To find the expecta- 
tion value of the total system energy E — {H), we use again 
Wick's theorem (see Appendix E) systematically. Explicitly, 
this energy functional is given as 

E{t) = Tr I [7^ + ^Uf^o, + Uj\ + [n+ Ufi., + Uj] f 

+ ^ mf'^)* + i [V„(., + Vs,] m*| = const. (58) 

That this is also a constant of motion follows straight from 
Eqs. (22,23) and make this a "derivable theory" according to 
Ref. [23]. For example, the same first order results can be 
found in Ref. [58], derived by a variational procedure. 



3. Entropy and individual occupation probabilities 

In a previous section Sec. II C, we have demonstrated with 
Eq. (18) that any admissible covariance matrix is necessarily 
of the following form 

G{t)^V{t)P{t)V{ty, (59) 
a3 = V{tya3V{t). (60) 

Furthermore, we have derived a non-linear equation of mo- 
tion (23) for this entity in Sec. II D. It is not self-evident that 
the aforementioned constraints on the structure of the density 
matrix are automatically preserved during the time evolution. 
Fortunately however, this is the case. This can be verified eas- 
ily by inserting Eq. (59) into Eq. (23). From this, one finds 



that all instantaneous eigenvalues of the density matrix are 
constants of motion 

P{t) = const., (61) 

provided that the instantaneous eigen basis evolves according 
to 

ifi^V{t) = T,V (62) 

In turn, this is compatible with Eq. (60) due to the symmetry 
of the self-energy in Eq. (36). 



ni. RESULTS 

A. Rescaling tlie kinetic equations to a prolate, quasi 
one-dimensional configuration 

General discussions on the properties of many-body 
physics are usually plagued by a very abstract language. In or- 
der to gain further insight into the complex non-linear physics, 
it is mostly necessary to examine spatially homogeneous sys- 
tems and assume stationary in the time-domain. This makes 
it possible to continue with analytical work and to capture 
the essential bulk physics of macroscopic systems. However, 
trapped quantum gases are different in many respects as the 
system size is of the same order as the coherence length, thus 
boundary effects are of equal importance. Moreover, it is dif- 
ficult to discuss succinctly the meaning of the thermodynamic 
limit in a quantum system with a mesoscopic particle number 
and a discrete energy spectrum. 




Figure 1: Schematic representation of the density distribution in a 
typical "cigar-shaped" trapped gas. The kinematic motion in the or- 
thogonal directions y,xis frozen out and presumably all the dynam- 
ics occurs along the weakly trapped s-axis. 



In the following sections, we will therefore examine nu- 
merically the properties of a quasi-one dimensional prolate 
system as shown in Fig. 1. While we have studied al- 
ready a three-dimensional isotropic configuration [71], pro- 
late "cigar-shaped" traps are currently at the focus of attention 
[40, 43, 72, 73]. Due to the reduction of the available phase- 
space volume, the role of quantum fluctuations becomes more 
pronounced in low dimensional systems. 

To be specific, we want to assume that the single particle 
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trap Hamiltonian of Eq. (4) is of the form 



H(x) 



2m 



A + \4o(x) - n, 



(63) 



Ho(x) = ^mu^x'^ + ^ma;_L^(2/2 + z^). (64) 

A measure of the anisotropy of this trap geometry is the as- 
pect ratio P = u)±/u. Ina very elongated, prolate configura- 
tion the perpendicular oscillation frequency u}± is much larger 
than the longitudinal circular frequency cu, thus [3^1. In 
order to reduce the kinetic equations also to a dimensionless 
form, we choose the ground state extension of the longitudinal 
harmonic oscillator gq = h/m, lu as the basic length scale, 
to = 27r/a; as the basic time-scale and Eq = w as the natural 
energy unit. 

The present formulation of the theory in terms of a generic 
set of basis states {|1)} is able to handle the dimensional re- 
duction from the general three-dimensional kinetic theory to 
the quasi one-dimensional situation quite easily, if we use an 
adapted trap basis in terms of longitudinal and perpendicular 
harmonic oscillator states |1) = \ji,ki)j_. By assum- 
ing that only the ground state components of the fields in the 
transverse directions are occupied, one can limit the evolution 
to one dimension effectively 

|a)=ai|l)0|O,O)^ (65) 
/=/i,2|l)(2|®|0,0)^(0,0|^, (66) 
m=mi,2|l)|2)«)|0,0)^|0,0)^. (67) 

In order to evaluate the matrix elements of the binary interac- 
tion potential further, we need the position representation of 
the normalized two-dimensional harmonic oscillator ground- 
state 



{y,z\0,0} = ipQ{y,z;t) = e 



TTttn 



"0 , (68) 



j dydz\ipo\'^ = 1, j dydz\ifo\' 



27ra; 



2 • 



(69) 



To compensate for the ground state energy of the two- 
dimensional harmonic oscillator, we are working here with 
an exphcitly time-dependent basis that removes that energy 
shift. We can now partially evaluate the matrix elements of 
the two-body interaction in the contact potential approxima- 
tion of Eq. (9). By assuming only ground states in the trans- 
verse direction, one obtains 



hujan — 
2 



6{x - y)n(x,t) and S(.x,y,t) ^ 5{x - y)T,{x,t). Fi- 
nally, this leads to the quasi one-dimensional self-consistent 
Hartree-Fock-Bogohubov (SCHFB) equations in the contact 
potential approximation 

idtx = 'n.{x,t)x{x,t), (71) 
idtG = i:{xut)G{xi,X2,t) - G{xuX2,t)'E{x2,t)\ (72) 

where local self-energies operators have been defined as fol- 
lows 



n(a;,t) 



ii^{x,t) n^(.T,t) 
-n*^{x,t) -u*^{x,t) J ' 



(73) 



+ g\a{x,t)f + 2gf{x,x,t) - /x, 



Il^{x,t) =gm{x,x,t), 
and 



j:{x,t) 



(74) 
(75) 



(76) 



^J^{x, t) = + 2g \a{x, t)\'' + 2g f{x, x, t) - m, 

(77) 

^A{x,t) = ga{x,t)^ + grh{x,x,t). (78) 

The numerical results that are presented in the following 
sections are based on typical *^Rb condensate data, i. e. , an 
atomic mass ms7 = 1.443110"^^ kg, axial and radial trap 
frequencies {u, vx) = (3, 800) Hz, which result in an aspect 
ratio P = 267, an axial harmonic oscillator ground state size 
ao = 6.2263 /xm and a 3d scattering length = 5.8209 rmi. 
This leads to an effective couphng constant g = 2pas/ao = 
0.4986. 



B. The Gross-Pitaevskii equation and the collective 
Bogoliubov excitations 

Before discussing the properties of the fully self-consistent 
equations, it is prudent to start with the most basic version of 
it. By disregarding / and m altogether, one obtains the Gross- 
Pitaevskii equation 



dx{l\x){2\x){x\S){x\4), (70) U^{x,t) 



(79) 



+ g\a{x,t)\''-fx, nGP = 0. (80) 



where we have introduced a dimensionless coupling con- 
stant g = 2Pas/ao. This corresponds to an effective one- 
dimensional interaction potential V^^\x) = hujaog S{x), 
which would produce the same matrix elements. In this hmit 
of a very locahzed, point-hke binary interaction, one finds fur- 
ther that the self-energy operator of Eqs. (24) and (31) be- 
come local operators in space and time , i.e., II{x,y,t) = 



From a stationary solution xo = (o^oi'^o)^' which is nor- 
malized to the total particle number N = J dx\aQ{x)\^, one 
obtains the chemical potential fi{N). In the interaction domi- 
nated mean-field regime, one may disregard the kinetic energy 
contribution altogether and one finds in the Thomas-Fermi ap- 
proximation fijp = (3.gA^/2)2/'^/2, the width of the conden- 
sate cctf = a/2/xtf and a healing length ^ = 1 /xtp. 
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1. Linear response analysis 

From a weak perturbation around the stationary solution 

a{x, t) — e^**''^[aQ(2;) + 5a{x, t)] or equivalently 

X{x,t) = e-'*'^^^^ [xo(x) + 5x{x,t)] , (81) 

one obtains the collective Bogoliubov excitations modes and 
linear response frequencies of Eq. (79). 

idtSx ^Y.''[x)5x{x,t)-5^lP, 



Elix) =- 



+ 2g|ao(x)|2 



(82) 

(83) 
(84) 



The energy spectrum for the linear response matrix must 
be gap-less (Eq = 0) as the GP-Eq. (79) is U{1) invaiiant 
under a global phase change of Oq. In other words, there ex- 
ists a degenerate manifold of ground state solutions and it re- 
quires no energy to transform one of them into another The 
infinitesimal generator of this phase rotation is the zero-mode 
P = (q!q, — Q!q)^ and it can be verified easily that S^P = 0. 

The initial perturbation Sx can also induce a small change 
in the particle number 5N — P^a^Sx- In turn, this leads to a 
small modifications of the chemical potential ^{N + 6N) = 
fi{N) + Sji at which the perturbed system is evolving globally. 
By taking this into account in the ansatz for the perturbation 
analysis in Eq. (81), one is able to cope with the secular terms 
that arise otherwise in a gap-less linear response analysis 



Sx{t) 



-it s' 



Sx{x,0) +itS^lP 



= iQ6N- iP 6^ + We-'' ^'ai"~^'^dx', (85) 

where we have used the completeness relation of Eq. (43) and 
defined Bogoliubov amplitudes as (5$ = Q^crsSx, as well as 
6x' = W^'a^Sx- The residual energy shift of the chemical 
potential is proportional to the inverse of the inertial mass of 
the collective ground state mode Sfi = 6N/M. Within the 
TF approximation, one finds for the mass parameter A/jf = 

(125iV)V3/g. 

In considering the similarities of the collective Bogoliubov 
excitations Eq. (82) and the quasi-particle modes of the gen- 
eralized self-energy matrix of Eq. (72), it is of utmost impor- 
tance not to confuse their different physical meaning. Quasi- 
particle modes address the questions of excitations of the 
quantum fluctuations above a static mean-field without con- 
sidering the back-action. Consequently, they do not have to 
be gap-less nor satisfy a Kohn theorem [74], but they are just 
a convenient basis to describe the quantum vacuum or the 
thermal excitations thereof. Thus, they should not be con- 
sidered as better or worse approximation of each other. Only 
the linear response analysis of the coupled system of Eqs. (7 1 , 
72) will be comparable to the collective Bogoliubov excita- 
tions of Eq. (82). It is straight forward to verify that the cou- 
pled system of Eqs. (71,72) are U{1) invariant, thus gap-less 
{Eq = 0). Moreover, the collective center of mass oscillation 
also decouples from the internal excitations and it evolves ex- 
actly with the harmonic oscillator frequency of i?i = 1. 
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Figure 2: The condensate density |a(a::)p as a function of the posi- 
tion X in units of the harmonic oscillator length for a particle number 
of A'^ = 10^. Only in the proximity of the Thomas-Fermi radius 
xtf = 9.07 the exact numerical solution (solid line) is distinguish- 
able from the Thomas-Fermi approximation (dashed doted line). The 
chemical potential energy in h.o. units fj, = 41.205 is also well ap- 
proximated by /iTF = 41.198. 




Figure 3: The normalized, even spatial components of the nega- 
tive energy Bogoliubov excitation mode W{E„ < 0) = [ue,«e]^ 
as a function of the position x in units of the harmonic oscilla- 
tor length. In particular, we show the three low energy modes 
2,4) = (0, —1.732, —3.167) and an arbitrarily chosen higher 
energetic mode Ego = —63.845 for comparison. 



2. Static mean-field and collective Bogoliubov excitations for 10 
particles 

In Figs. 3 and 4, we depict a few selected collective excita- 
tion modes of the GP-equation corresponding to the static so- 
lution with N — 10'^ particles shown in Fig. 2. The reflexion 
symmetry of the harmonic trap is carried over to the mean- 
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Figure 4: The normalized, odd spatial components of the negative 
energy Bogoliubov excitation mode W{En < 0) — [uo,Vo]^ as 
a function of the position x. Depicted are the three low energy 
modes _B„^(i 3 5) = ( — 1, —2.451, —3.882) and an arbitrarily cho- 
sen higher energetic mode Esi — —64.723 for comparison. 



field state. Consequently, one can also classify the collective 
Bogoliubov modes according to even and odd parity modes. 
Conventionally, one introduces also hole u{E < 0) and parti- 
cle v{E < 0) amplitudes as components of W{x, En < 0) = 
[u{x, E), v*{x, E)]^ . The particle character of v{E < 0) be- 
comes visible in the higher energy excitation modes i^so.si 
where only the particle-like amplitudes reach far outside the 
spatial extend of the condensate and oscillate between the 
classical turning points of fictitious particles with energy \E\. 
In contrast, hole-like excitations always remain localized on 
the site of the condensate wave-function and their amplitudes 
decrease with decreasing excitation energy. 



3. Chemical potential and collective Bogoliubov excitation 
energies for a full range of particle numbers 

The general character of the solutions depends parametri- 
cally on the particle number. In Fig. 5, we have plotted the 
chemical potential /i(A^) as well as the collective Bogoliubov 
energies E{N) as a function of the particle number N . The 
double logarithmic representation of the chemical potential 
shows clearly a linear slope, hence reproduces the power law 
dependence that is found within the Thomas-Fermi approxi- 
mation fiTF{N) = (3giV/2)2/3/2. The collective excitation 
frequencies are identical to the harmonic oscillator spectrum 
for small particle numbers. However with increasing particle 
number, the energy density of states visibly increases as the 
level spacing is reduced. By focusing on the two lowest en- 
ergy excitations, i. e. , the Goldstone mode Eq = and Kohn 
mode El = 1, one finds that those fundamental symmetry 
generating modes are unaffected by the interactions. 




n t t t t t t — \ 

10"^ 10° 10^ 10^ 10^ 10* 10^ 
N 



Figure 5: The upper part shows the chemical potential n{N) of the 
GP equation as a function of particle number (solid line). For com- 
parison, we present also the Thomas Fermi approximation fiTp{N) 
(dashed line) and the harmonic oscillator ground state energy of 1 /2 
(dashes doted line). In the lower part, the collective Bogoliubov en- 
ergies E are shown versus particle number. All units are given in h.o. 
energies. 



C. Ground state solutions of the self-consistent 
Hartree-Fock-Bogoliubov equations 

In this section, we will discuss the self-consistent solution 
of the stationary Eqs. (71) and (72), i. e. , 

G = nx, = SG-GS^ (86) 

In particular, we are interested in the lowest energy configu- 
ration of the system. Thus, no quasi-particle modes shall be 
occupied and according to Eq. (21), we construct the quantum 
vacuum only from the negative energy states W- — We<o 
with dim[Ty_]=2n x n 

G = W_Wl. (87) 

With an iterative procedure, one can solve for the mean- 
field amplitude q;q assuming static / and m. In a second 
step, one has to find the quasi-particle modes and construct 
a density matrix from it. This procedure has to be contin- 
ued until convergence is reached. We have computed the 
self-consistent solutions for a range of particle numbers N = 
(10°, 10\ 10^, 10^, ID", 10^). For a particle number N < 1, 
this extended mean-field theory is neither physically meaning- 
ful, nor did it lead to converging solutions any longer. 



1. Static mean-field and quasi-particle modes for 10^ particles 

For a particle number N = 10^, we show in Fig. 6 
the mean-field density as a function of position. Due to a 



11 



|a(x) 



80 
70 
60 
50 

2 

40 
30 
20 
10 




































































































































































-15 -10 -5 5 1 


15 



1 

0.5 
u'(x) 
-0.5 
-1 



1 1 1 


1 1 1 


1 1 1 






















»i;WWWWVW\ 










































-15 -10 -; 


) 5 10 15 




Figure 6: The self-consistent mean-field density (solid line) 

and Thomas-Fermi approximation (dashed doted line) as a function 
of the position x in units of the h.o. length. 
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Figure 7: The normalized, spatial components of the negative en- 
ergy Bogoliubov excitation mode W{En < 0) = [ur,v*.]^ as a 
function of the position x in units of the harmonic oscillator length. 
In particular, we show the three low energy modes i?„={o,2.4) ~ 
(—3.356, —6.421, —8.636) and an arbitrarily chosen higher ener- 
getic mode Eso = —64.819 for comparison. The low energy modes 
are localized on the right side of the condensate and almost degener- 
ate to the modes localized on the other side depicted in Fig. 8. 



repartitioning of particles between condensate and the non- 
condensed fraction, there are now fewer particles in the mean- 
field component. Thus, the difference to the Thomas-Fermi 
approximation is more visible than in Fig. 2. For a total parti- 
cle number of TV = 10^ = + N, we find a {N'^"^ , N) = 
(984.83, 15.17) and chemical potential of ^ = 40.61 in h.o. 
energy units. 

The quasi-particle modes are depicted in Figs. 7 and 8. 
While the collective Bogoliubov modes of Figs. 3 and 4 can be 



Figure 8: The normalized, spatial components of the negative en- 
ergy Bogoliubov excitation mode W{En < 0) = [ui,vi]^ as a 
function of the position x. Depicted are the three low energy modes 
_B„^(i 3 5) — (—3.3556, —6.421, —8.636) and an arbitrarily cho- 
sen higher energetic mode Esi — —65.689 for comparison. The 
low energy modes are localized on the left side of the condensate. 



characterized with a definite even or odd parity, this is seem- 
ingly not the case here. The low energy quasi-particle modes 
are localized on left and right sides of the condensate and have 
more of the character of the single particle Hartree-Fock ex- 
citations discussed in Eqs. (32) and (49), where the potential 
energy Vho + 2Z^y(c) — /i has a double minimum at ixxp. How- 
ever, one has to take into account also that the low lying modes 
are energetically degenerate and that the higher energy exci- 
tations do exhibit a definite parity. Thus, one could construct 
quasi-particle modes with definite parity by symmetrizing or 
anti-symmetrizing them. We have deliberately chosen not do 
so in order to permit the occurrence of a reflexion symmetry 
breaking. For example, spatial symmetry breaking is observed 
in deformed nuclei [58] or the mixing of two component Bose 
gases [75]. Nevertheless, we did not find such a symmetry 
breaking behavior here. 



2. Chemical potential and quasi-particle energies for a full range 
of particle numbers 

In Fig. 9, we present flie chemical potential ii{N) as well as 
the quasi-particle energies E as function of the total particle 
number N . One can see very clearly the energy gap Eq ^ Q 
and the absence of the Kohn mode Ei ^ 1, which should 
occurs only for the collective excitations of system as whole. 
Moreover, one can observe how the low lying quasiparticle 
modes begin to coalesce with gradually increasing particle 
number. 
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Figure 9: Upper part: the chemical potential fi{N) (solid line) of 
the SCHFB equations, the TP approximation /itf (dashed line) and 
the h.o. ground state energy 1 /2 (dashed dotted line) versus particle 
number A'^. Lower part: quasi-particle energy spectrum E(N) versus 
particle number A'^. 




Figure 10: Coherent contribution of the single particle density matrix 
Z''^' (xi, X2) as a function of the spatial coordinates for a particle 
number of N = 10'^. The coherence extends in the diagonal as 
well as in off-diagonal directions up to the TP radius. For a real 
valued mean-field q(x), the coherent contribution of the pairing field 
m''^' (xi,X2) is also represented by this figure. 



D. The single particle density matrix and the pairing field 

The mean-field amplitude a{x), the total single particle 
density f{xi,X2) = f'''^H^iT^2) + fixi,X2), as well as the 
pairing field 'm{xi, X2) have been the central concepts of the 
present analysis. Thus, we will present in the following pic- 
tures instances of their spatial representation for a particle 
number of iV = 10'^. This gives a good qualitative impression 




Figure 11: Non-condensate density matrix f(xi,X2) as a function of 
the spatial coordinates for a particle number of A'^ = 10"^. 



of the universal features of the ground state. Results obtained 
for different particle numbers are similar in appearance and 
we will discuss the quantitative differences next. 



1. Spatial representation for 10'^ particles 

In Fig. 10, we show the coherent contribution to single par- 
ticle density matrix 2:2) = a*{x2)a{xi). As it is 
constructed from the order parameter a{x), it has the full 
ODLRO [20, 21], which extends over the complete system 
size. For the one-dimensional trap that we consider in here, 
the Hamiltonian operator is real valued, thus the ground state 
solution of the mean-field a{x) is a purely real quantity, too. 
Consequently, m'^'^^(a;i, X2) = a{x2)a{xi) is identical to 
f'^''\xi,X2) as depicted in Fig. 10. 

In contrast to the full ODLRO which is found in the mean- 
field component, one finds that the non-condensate density 
/(xi, X2) is predominantly localized along the coordinate di- 
agonal and exhibits spatial variation only due to the external 
confinement with a trapping potential. In the off-diagonal di- 
rection this order parameter has only a very short range which 
is determined by the binary interaction. This is depicted in 
Fig. 11. 

Interestingly, one finds also for the pairing field m{xi, X2) 
a very similar spatial behavior. However, while a single par- 
ticle interpretation is sufficient to understand the behavior of 
the normal density matrix, it is necessary to use two-particle 
physics in the pairing field of Fig. 12. The strong negative cor- 
relation along the diagonal show that there is a reduced like- 
lihood of finding two particles at the same location. Again, 
this likelihood is modulated by the density of particles in trap. 
However, the degree of anti-correlation drops off quickly with 
increasing distance between the positions at which particles 
are extracted. 
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Figure 12: The negative pairing field —m{xi,X2) as a function of 
the spatial coordinates for a particle number of A'^ = 10'^. 



2. Diagonal and off-diagonal elements of f and m for a full range 
of particle numbers 

In this section, we have compiled the quantitative results for 
the spatial variation of the diagonal /(x, x), ■m{x, x) and the 
off-diagonal f{x, —x), fh{x, —x) elements of the normal den- 
sity matrix and the pairing field, respectively, for a full range 
of particle numbers N — (lO'^, . . . , 10^). The double loga- 
rithmic representation used in Figs. 13 and 14 reveals clearly 
that there is a separation of the bulk physics in the center of 
the trap and the physics dominated by the boundary at the rim 
of the condensate. We have verified this separation of scales 
by turning of the trap potential. In this case, one recovers 
the homogeneous limit. In order to make a quantitative com- 
parison with the trapped system, we have chosen a homoge- 
neous mean-field density n = N/2L = 82.62 that matches 
the mean-field density at the center of the trapped gas for a 
particle number of = 10"^. This leads to very similar chem- 
ical potentials for the homogeneous system of /ihom — 40.47 
and for the trapped gas of fiaap — 40.61, respectively. The 
dashed dotted lines in Figs. 13 and 14, do represent the ho- 
mogeneous results and compare very well with the trapped 
gas for N = 10"^. This proves that the local density approxi- 
mation (LDA) yields a good approximation for the transverse 
coiTelation length. A detailed comparison of the critical expo- 
nents and their dependence on temperature is cuiTently under 
investigation. 



E. Effective coupling constants 

In Sec. II F, we have discussed formally the renormaliza- 
tion of the binary interaction to an effective T-matrix, which 
takes place in any self-consistent calculation. In particular, we 
upgraded the anomalous potential of Eq. (53) and related it 
to an anomalous many-body T-matrix. As we have calculated 




0.1 0.2 0.4 



Figure 13: Diagonal f{x,x) and off-diagonal \ f{x, -~x)\ elements 
of the normal density matrix versus distance x for a full range of 
particle numbers TV — (lO", . . . , 10^). The individual results can 
be identified easily by the spatial extension that grows proportional 
with the particle number. The dashed-dotted line shows the result for 
a homogeneous gas corresponding to a trapped gas with A'^ — 10^ 
particles. 
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Figure 14: Diagonal —m{x,x) and off-diagonal \m{x, —x)\ ele- 
ments of the pairing field versus distance x as in Fig. 13. 



the SCHFB equations in here, we can now revisit this question 
and find out what the effective renormaUzed coupling constant 
is 

TA{2n,x)m^''\x,x) = g[m^''\x,x) +m{x,x)], (88) 
T^i2^i,x)f'^''\x,x)^2g[f<^''\x,x)+f{x,x)]. (89) 

From the numerical results of Fig. 15, we can draw three con- 
clusions: a) the effective coupling is position (or momentum) 
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Figure 15: The effective coupling constants TAr(2/i, x) and 
rA(2/i, x) as a function of position. Each subplot contains six in- 
dividual lines corresponding to the full range of particle numbers 
— (lO", . . . , 10^) and a flat line, which the represents the bare 
coupling constant g = 0.4986. The individual results can be identi- 
fied easily by their spatial extension, which increases proportional to 
the TF radius or particle number, respectively. 



dependent and decreases rapidly outside of the range of the 
condensate; b) the effective coupling constant is in general 
less than the bare interaction constant, which may be under- 
stood in terms of a second order perturbation theory; c) the 
effective coupling constant at the center of the trap becomes 
gradually less for smaller particle numbers. 



F. First and second order correlations functions 

1. Spatial representation for 10^ particles 

Quantum fluctuations around the classical mean-field am- 
plitude are the central topic of this article. In the previ- 
ous sections, we have examined the single particle density 
fx = Z''^' + /, which is directly an observable quantity and the 
pair correlation function iri{xi, X2), which is not. The ques- 
tions of how to quantify and to measure quantum correlations 
has always been a central theme for any quantized field the- 
ory, whether in condensed matter physics [20, 21, 32] or in 
quantum optics [76, 77]. 

In essence, first order coherence is measured by the corre- 
lation function 



(90) 



y^n{xi) n{x2) 
_ .f^''Hxi,X2) + f{xi,X2) 
\/n{xi) n{x2) 

where n{x) = f^'^^ (x, x) + /(x, x) denotes the total density. 
Primarily, it is sensitive only to spatial phase correlations or. 
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Figure 16: First order correlation function g^^^ {x\ , X2) versus spatial 
position for 10'^ particles. 



in other words, off-diagonal order. This can bee seen eas- 
ily by disregarding the quantum depletion / for the moment. 
In the case of a mean-field a{x) = |Q;(a::)|e^**'^^^, it is only 
proportional to the phase gradients <i>(x) — $(0) + 5^{x). 
Moreover, if a{x) is a static ground state without a phase gra- 
dient (irrotational), one finds the definition of full coherence, 
\.&.,g'^^\xi,X2) = 1. 

In Fig. 16, we have evaluated the first order correlation 
function for a trapped gas with N = 10'^ particles. By def- 
inition, it must be exactly 1 along the diagonal. As expected, 
there is only a very weak influence of the quantum fluctua- 
tions noticeable, since f'^'^\xi,X2) S> f{xi,X2), which can 
be seen explicitly in Fig. 6 and Fig. 11. Consequently, one 
finds that first order coherence is not a very sensitive probe for 
the quantum aspects of a field. For example, sending a clas- 
sical optical field through an semi-opaque or noisy medium 
immediately leads to a reduction coherence, which is as such 
a purely classical phenomenon. 

Actually, a sensitive probe for the quantum nature of a field 
is the second order coiTelation function g'^'^\ If g^^^ basically 
responds to uncertainty in the phase quadrature, then g^^^ is 
affected by density fluctuations aj,aj. — n{x) + 5n{x). In- 
tuitively speaking, this is the conjugate variable to the phase 
gradient. More succinctly speaking, this can be examined by 
studying a Heisenberg uncertainty product (see App. A) for 
all the field quadratures as in Eq. (46). Explicitly, it is defined 
as the normal ordered density-density correlation function: 



9^ '{XI,X2) 



-A / 

{xi)n{x2) y 



n{xi)n{x2) n{xi)n{x2) 
f^''\xi,X2)* f{x2,xi) + m^''\xi,X2)*m{x2,xi) 

+ f{xi,X2)f{x2,xi) + m{xi,X2)*m{x2,xi)^, (91) 



which is shown in Fig. 17 for 10'^ particles. While g^^-* it is 
mostly equal to 1, it can be seen clearly that along the diagonal 
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Figure 17: Second order correlation function g^-^\xi, X2) versus 
spatial position for 10'^ particles. 



g^^-* {x, x) < 1. This is a unique signature a non-classical state 
of the quantum field with a sub-poisonian statistic of number 
fluctuations or, in other words, a number squeezed state [34] . 
The suppression of number fluctuations has its origin the bi- 
nary interaction potential of Eq. (1). By reducing the local 
density fluctuations, one can remove more energy than what 
is regained by the kinetic energy due to the increased phase 
uncertainty. At its extreme this is also the physical mecha- 
nism that leads to the Mott phase transition [11-13, 78] where 
on-site interactions compete with nearest neighbor tunneling. 



2. Diagonal and off-diagonal elements of the first and second 
order correlation function for a full range of particle numbers 



Finally, we summarize the results for the first and second 
order correlation function in Figs. 18 and 19 for the full range 
of particle numbers N = (10", . . . , 10^). As before, we have 
also computed the result for a homogeneous gas in the local 
density approximation, which corresponds to a trapped gas of 
10'^ particles. The insets in the pictures magnify both of those 
curves in the central region where they do agree very well. 
However, at the rim of the condensate the trapped gas does 
exhibit features that are absent in the LDA. 

It is physically most important to see in Figs. 19 that the 
suppression density fluctuations becomes stronger for smaller 
particle numbers. At its limit, this anti-correlation of tow 
bosonic particles leads to an effective fermionization and is 
the hall mark of the Tonks-Girardeau regime [18, 37, 50]. 
If this argument is applied to the situation of an inhomoge- 
neous trapped gas, then this means that the number squeezing 
is larger at the rim than in the center of the cloud. 
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Figure 18: Off-diagonal elements of the first order correlation func- 
tion g^^\x, —x) versus spatial position. Each subplot contains six 
individual lines corresponding to the full range of particle numbers 
N = (10", . . . , 10^) and a dashed dotted line, which represents the 
homogeneous gas result. The individual results can be identified eas- 
ily by their spatial extension, which increases proportional to the par- 
ticle number. 
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Figure 19: Diagonal {x,x) and off-diagonal g^ '{x,—x) ele- 
ments of the second order correlation function versus spatial posi- 
tion. Parameters and legend as in Fig. (18). 



IV. CONCLUSIONS AND OUTLOOK 

In this article, we have given a basic overview of the 
premises and concepts of a number-symmetry broken non- 
equilibrium kinetic theory of a trapped bosonic gas. By 
extending the successful mean-field concept of the Gross- 
Pitaevskii equation with the effects of non-local, two particle 
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quantum correlations, one obtains a renormalized binary in- 
teraction T-matrix and allows for the dynamic establishment 
of non-classical many-particle quantum correlations. At very 
low temperatures, we can disregard the equilibrating effects 
of elastic two-particle collisions, in contrast to previous work 
[27, 28, 71]. In this limit, we have proven that the important 
physical constants of motion, hke particle number, energy or 
the entropy are conserved. Obviously, the inclusion of colli- 
sion processes is desirable from a fundamental point of view, 
as this will break the micro-reversibihty and lead to a thermal 
relaxation towards the most probable distribution as dictated 
by thermodynamics. In preliminary study [31], we have al- 
ready studied the effects of collisions and the consequences 
of various approximations in the absence of a mean-field in a 
model system, but further work is necessary. 

In the second section of the paper, we have focused on the 
specific properties of the zero-temperature ground state corre- 
lations of a trapped, quasi-one dimensional bosonic gas. This 
maximizes the role played by quantum fluctuations, due to 
the dimensional reduction of the available phase-space vol- 
ume. With a fully self-consistent numerical calculation, we 
have evaluated, the mean-field ampUtudes, the quantum de- 
pletion, the pairing field as well as the first and second or- 
der ground state quantum correlations for a full range of par- 
ticle numbers N = (10°, . . . , 10^). Most interestingly, we 
do get a strong suppression of the density fluctuations or, in 
other words, an enhanced number squeezing with decreasing 
particle density. This generic feature is in general agreement 
with the predictions that are found with exactly solvable one- 
dimensional models such as the Tonks-Girardeau gas or the 
Bose-Hubbard gas on an optical lattice. A detailed analytical 
comparison with those models [37-39] and the experimen- 
tal results [40, 43] is currently work in progress. In conclu- 
sion, we find that the general non-equilibrium kinetic theory 
also reproduces the ground state correlations of a quasi one- 
dimensional bosonic gas well. 
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Appendix A: CAUCHY-SCHWARTZ INEQUALITY 

For a positive semi-definite density operator a and an arbi- 
trary operator L it follows that the expectation value 

(LL^) = TT{aLL''} > (Al) 

is never negative. Consequently, the co-variance matrix G of 
Eq. (14) must be positive semi-definite u^Gu > 0, as well. 
This can be easily seen, by considering a linear combination 
of two arbitrary operators A and B, i.e. , L = aA + (313. By 
minimizing the positive expression Eq. (Al), one obtains the 



Cauchy-Schwartz inequaUty as 

(lit ) {BB'i ) > (Bit ) (ifit ) . (A2) 

In particular, for the special choice of A = 5di and B = da\, 
this implies that the magnitude of the anomalous fluctuations 
is Umited by 

(l + /ii)/22 > |mi2p. (A3) 



Appendix B: PAULI MATRICES 

In this article, we use the foUowing standard representation 
for the PauU matrices: 

f t\ . f -t\ f t \ 

(Bl) 

They satisfy the ordinary commutation relation of an angular 
momentum operator [ci , 0-2] = «f 3 and all cyclic permuta- 
tions thereof. If n is the dimension of the vector space, then 
cTft acts in a 2n-dimensional symplectic vector space [58]. 



Appendix C: CANONICAL TRANSFORMATIONS 

A canonical transformation is an inhomogeneous linear 
combination of creation and destruction operators that pre- 
serves the commutation relation [58]. In particular, if a and 
flt denotes a pair of hermitian conjugated bosonic operators, 
such that 

[a„al] = 5i,2, (CI) 

then any affine linear transformation defines a new set of op- 
erators b and b by 

In an n-dimensional vector space, T represents a 2 n x 2 n 
dimensional matrix and is a 2 n dimensional vector. Such 
a transformation is canonical if the new pair of operators also 
satisfies the commutation relation: 

[6i,H=<5i,2. (C3) 

More specifically, the transformation is unitary canonical if 
the new operators are hermitian conjugate pairs, i. e. , 6 = 6^. 
By inserting Eq. (C2) into Eq. (C3), one finds that the transfor- 
mation matrices are a representation of the symplectic group 

Sp{2n): 

Ta3T^ = a3. (C4) 

In addition, it can be shown that T* = cti Tcti and = 
as TVs. 
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Appendix D: QUANTUM LIMIT FOR THE GROUND STATE 
CORRELATIONS 

The state of the interacting many-body system is described 
within the set of approximations by a mean-field amplitude a 
and a generalized G matrix. At finite temperature or out of 
equilibrium, there is a finite occupation number of particles of 
positive energy excited states. However, in the lowest energy 
configuration all allocatable real particles occupy the mean- 
field amplitude and the generalized density matrix only holds 
the vacuum. Then, the following idem-potency relation holds 
for the density matrix 

Ga3G + G = 0. (Dl) 

This follows straight from Eq. (21) when P+ = 0. Conse- 
quently, this gives a restriction for the components of the den- 
sity matrix /, m. In particular, one finds that the vacuum de- 
pletion is completely determined by the pairing field 

/=^(^l+4mmt-l), (D2) 

fm = mf*. (D3) 

It is important to note that this relation also holds out-of- 
equilibrium as long as thermaUzing colUsions can be disre- 
garded. 

Appendix E: A GENERALIZED WICK'S THEOREM 

GauBian fluctuations around a well defined mean value are 
the key assertion to apply Wick's theorem [56]. This is a set 
of rules to efficiently evaluate quantum averages for multiple 



operator products as 

{i^l^2-. 4l){c.,c.' J, fh} (El) 

In this average, for example, the operator tjji represents either 
an operator or a|. 

First, the displacement rule shifts any operator ifji by its c- 
number expectation value ipi which is either or al, and re- 
places the quantum average by an average that has zero mean 
values: 

= ((^1 + ^/'l)(lA2 + ■ • ■ ii^n + V'i)>{0,0./,m}- 

Second, after expanding the multiple products, one can dis- 
card all averages that involve an odd numbers of operators: 

(V'lV'2..-V'26+l)(o,0,/,m} = 0- (E3) 

And third, for the remaining averages, one can use the 
GauBian factorization rule: 

(V'lV'2...V'2.){o,0,/,m} = (E4) 
= (V'lV'2)|o,0,/»(^3 • • • ^2«){o,0,/,m} + 
+ (V'lV'3){0,0,/>}(^2V'4 . • . V'2.){o,0,/,m} + 

+ (V'lV'2s){0,0,/,ifi}(V'2 • ■■i'2s-l)[o,0,f,fh}- 

By proceeding recursively, one has finally evaluated the com- 
plete multiple operator average Eq. (El). 
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